Deconfinement transition dynamics and early thermalization in QGP 
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We perform SU(3) Lattice Gauge Theory simulations of the deconfinement transition attempt- 
ing to mimic conditions encountered in heavy ion collisions. Specifically, we perform a sudden 
temperature quench across the deconfinement temperature, and follow the response of the system 
in successive simulation sweeps under spatial lattice expansion and temperature fall-off. fn mea- 
surements of the Polyakov loop and structure functions a robust strong signal of global instability 
response is observed through the exponential growth of low momentum modes. Development of 
these long range modes isotropizes the system which reaches thermalization shortly afterwards, and 
enters a stage of quasi-equilibrium expansion and cooling till its return to the confinement phase. 
The time scale characterizing full growth of the long range modes is largely unaffected by the con- 
ditions of spatial expansion and temperature variation in the system, and is much shorter than the 
scale set by the interval to return to the confinement phase. The wide separation of these two scales 
is such that it naturally results in isotropization times well inside 1 fm/c. 

PACS numbers: 11.15.Ha, 12.38.Gc, 25.75.-q, 25.75.Nq 



I. INTRODUCTION 

Recent heavy ion collision experiments at RHIC have 
produced a wealth of data on hadron spectra and their 
anisotropics, in particular the magnitudes of radial and 
elliptical flows. This data reveals, perhaps unexpect- 
edly, coherence in particle production and strong collec- 
tive flow phenomena. It turns out that about 99% of 
the single hadron data {px < 1.5 GeV) are very well 
described by the hydrodynamics of a near-perfect fluid, 
provided the initial condition of very rapid thermaliza- 
tion (in ^ 0.5/to/c) is introduced 0] - 0]. This strongly 
indicates that the quark-gluon plasma (QGP) formed at 
RHIC energies is a strongly coupled fluid. 

At asymptotically high temperatures above the decon- 
flnement Tc, where the running coupling g{T) is small, 
QCD is well-described as a gas of weakly coupled quasi- 
particles, and has been much studied by perturbative 
techniques. At the energy densities achieved in the high 
energy heavy ion collisions, however, perturbative treat- 
ment of the equilibration process appears not to be ap- 
plicable. Various estimates of thermalization times based 
on perturbative scattering processes have been obtained, 
e.g. in the so called parton-cascade approach to the time 
evolution of hard partons, or the bottom- up scenario jj, 
0, 01 incorporating saturation picture (see 0, |1| for 
review) initial conditions. They all result into thermal- 
ization times much longer than those needed by the hy- 
drodynamical simulations. It has been argued that this is 
a generic feature of any dynamical evolution based only 
on perturbative scattering processes 

Even within a weak coupling analysis, however, non- 
perturbative effects may contribute to the dynamics. 
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It has been pointed out that, in a plasma with an 
anisotropic hard parton distribution, instabilities may 
develop in soft gauge field modes generated within the 
linear response approximation ^3 " been ar- 

gued that these small-amplitude unstable modes are not 
stabilized by (non-Abelian) non-linearities; and thus can 
grow to amplitudes large enough to contribute 0(1) frac- 
tion to the total energy density, and drive isotropization 
of the hard modes faster than any hard collision equili- 
bration time scale . 

Various investigations of the evolution of these semi- 
classical instabilities have been carried out probing be- 
yond the linear regime, a nd, in particular, within the 
hard- loop effective action [13 - |23|- They generally in- 
dicate that such instabilities indeed persist when non- 
Abelian non-linearities are taken into account within the 
weak coupling regime. At some point, however, other 
non-perturbative dynamics at strong effective couplings 
defined at scales appropriate to the nonlinear interaction 
of such growing long range modes must enter. Still, con- 
sideration of such potential instabilities properly point to 
a basic underlying question which should be addressed 
from a more general point of view. 

When systems are driven far from equilibrium by sud- 
den changes in external conditions, the approach to a new 
equilibrium state involves, in general, complex nonequi- 
librium processes. Such situations arise, for example, 
in a 'quench' to the metastable region across a first- 
order transition boundary, or from a one-phase region 
to the multiphase coexistence region across a second- 
order transition. The study of the dynamics of such 
far-from-equilibrium processes is still at an early stage 
of development. Nonetheless, from many studies, mostly 
in condensed matter physics systems, two broad classes 
of responses have been roughly identified. 

Immediately after a rapid quench, the state of the sys- 
tem, as characterized by the appropriate order parame- 
ter, is nearly identical to the state before. The fields must 
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then find some way to adjust toward values appropriate 
to the conditions after the quench. One way the system 
may decay towards equihbrium is by the excitation of fi- 
nite amplitude localized fluctuations that may grow or 
coalesce as in nucleation, or interact in other ways over 
time. This typically indicates that the system finds it- 
self in some sort of metastable state. Another way is by 
the immediate development of a spatially modulated or- 
der parameter whose amplitude grows continuously from 
zero throughout the sample. Spinodal decomposition is 
a prime example of this type of response, and the term 
is often used loosely to generically denote such globally 
unstable behavior. It should be pointed out that the 
boundary between these two rough classes is not sharp 
in all cases. 

At RHIC heavy-ion collisions in the central rapidity 
region achieve the sudden deposition of energy densi- 
ties reaching ~ 30GeV/ fm^. The first most basic ques- 
tion that must be posed then is: which general type of 
dynamic response is characteristic of a rapid transition 
across the confinement- deconfinement boundary in QCD? 
In the case of heavy- ion collisions, the question must be 
further qualified by the inclusion of the effects of rapid 
expansion and temperature variation, as well as finite 
volume. This is the issue we explore in this paper. 

To get some intuition, we first investigate this question 
briefly within an effective action approach in Section 2. 
Indeed, much of the current understanding of the early 
time evolution of systems out of equilibrium has been 
obtained by investigating classes of stochastic equations 
that are natural dynamical (time dependent) generaliza- 
tions of the Landau-Ginzburg (LG) effective action mod- 
els of the static (equilibrium) theory (H^, In the 
case of the confinement - deconfinement transition, the 
relevant order parameter is the Polyakov loop (Wilson 
line), and LG effective actions for it have been consid- 
ered in |2^ - [25 . Corresponding dynamical model gen- 
eralizations can then be used to examine the question 
posed above. We consider the predictions of such a model 
briefly in Section 2 below. 

Though the effective model approach often proves valu- 
able, what is ultimately needed is an ab initio treatment 
in the full nonperturbative formulation of the exact the- 
ory, i.e. lattice gauge theory. Unfortunately, there is no 
established simulation formalism for directly extracting 
physical properties in non-equilibrium real-time evolu- 
tion in quantum field theory. What one can do, how- 
ever, is mimic Minkowski real-time dynamics by Glauber 
stochastic dynamics evolution. Thus starting with the 
system in thermalized equilibrium, one performs a tem- 
perature quench and follows the system, over successive 
simulation sweeps, on its path toward regaining equilib- 
rium. Though this cannot be directly identified with the 
exact real-time evolution, it is known from many studies, 
mostly of condensed matter systems, to accurately reflect 
it. At the very least, the method provides a consistent 
picture of the basic features of the system's real-time re- 
sponse. It has been extensively and successfully used for 



many systems exhibiting, in particular, first-order tran- 
sitions. Indeed, in studies of binary alloys and binary 
fluids it has been found to reproduce the experimentally 
observed behavior in quantitative detail '2^ . In the case 
of gauge theories, the method was first used in the pio- 
neering studies in 113, lll,!!^. More, recently, such stud- 
ies of the dynamics of phase transitions in spin and gauge 
theory systems were undertaken and much extended in 

In this paper, building on these previous gauge theory 
studies, we consider SU{3) gauge theory under condi- 
tions mimicking the situation encountered in heavy-ion 
collisions. This we do in Section 3 which constitutes 
the main part of the paper. Specifically, we examine 
the response after a sudden quench into the deconfine- 
ment phase under varying conditions of spatial expan- 
sion and temperature variation. In this first investiga- 
tion, we consider only pure SU{3) gauge theory, i.e. no 
quarks, as the deconfinement transition is driven by glu- 
onic dynamics. In simulation measurements of structure 
functions and Polyakov loops, we follow the evolution of 
the system from the quench till its return to the con- 
finement phase. The main outcome of our study is that 
two relevant, widely separated time scales emerge. First, 
a strong and robust signal of rapid growth of very long 
range modes is observed after the quench. Most impor- 
tantly, the time scale of full development of these low 
momentum modes is essentially unaffected by the pres- 
ence, over a wide range of parameters, of spatial expan- 
sion and temperature variation in the system. The de- 
velopment of these modes 'isotropizes' the system, which 
reaches full quasi-equilibrium shortly afterwards as sig- 
naled by the full decay of the structure function. It thus 
enters a stage of expansion and cooling, characterized by 
a second, much longer time scale, till its return to the 
confinement phase. The wide separation of these two 
scales accords well with what is seen in heavy-ion colli- 
sion experiments. Having arrived at this robust qualita- 
tive picture, we also attempt to make a more quantita- 
tive comparison of what is seen in these simulations to 
hydrodynamical phenomenology (subsection 3.3). There 
are certainly uncertainties here, not least of which is the 
fact that we do not include fermions which can make a 
significant contribution especially at the late stage near 
the return to confinement. Still, one finds that the sep- 
aration of scales is such that, for any reasonable choice 
of parameters, isotropization times well inside 1 fm/c 
result naturally. 

Finally, in Section 4 we briefly discuss our conclusions 
and directions for further work. 



II. EFFECTIVE ACTION MODELS 

Effective action models allow one to build simple the- 
ories of initial time evolution of systems driven out of 
equilibrium by a quench. Such theories can be build for 
general classes of models. A simple linear theory (Cahn- 
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Hilliard) was first proposed for models of type B with 
conserved order parameter [3^. The generalization to 
models of type A with non-conserved order parameters 
is straightforward, see , fo r example, [s^,!!^- For a gen- 
eral overview see e.g. |2lL |22 |. 

Though such effective models are not our primary fo- 
cus in this paper, it is useful to briefly consider them as 
they provide useful insight into the possible behavior of 
the system that can be checked against the outcome of 
simulations in the actual gauge theory. For SU{3) gauge 
theory the low energy degrees of freedom are represented 
by Polyakov loops, and standard potential models for 
them are known and well studied |U |3 HEl • Adopting 
the potential for the Polyakov loop I (a complex quantity) 



v(0 
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the coupled set of Langcvin equations is 

T], and its c.c. 



m _ 5S_ 

dt^^ 51* 



(1) 
(2) 



the homogeneous part of these equations. Its eigenvalues 
are 



Wl,2(fc) = -r(ifc2-^Ci±|c2|), 



and the eigenvectors are 

<5'l,2 = 

The solution is 



±C2/|C2| 
1 



V v{k,t) ) 



(7) 



(8) 



(9) 



Here, an analysis of modes similar to that of the real 
case (spin systems, S'C/(2)) (js^ - [13) can be applied. If 
ci ± |c2| < 0, one will observe exponential growth. Spin- 
odal decomposition-like behavior then results if we have 
exponential growth in at least one exponent, i.e. when 
ci < \ci\. The structure function (connected Polyakov 
2-point function) 



Here S is the standard Landau-Ginzburg action 



S{k,t) = ( u{k,t)v{-k,t) 



(10) 



S 



(3) 



r is the response coefficient, which defines the relaxation 
time scale of the system, and 77 is a noise term. We will 
ignore the noise term, since it can be shown that it does 
not affect the resulting rate of growth of fluctuations ■ 
We are interested in the fluctuations of the Polyakov 
loop I around some average value 1^: 



l{f,t) = lo + u{r,t), 



(4) 



Using and the system of equations for the fluc- 
tuations is 



— = -r 
dt 



-iv^U + (ciM -I- C2U* + C3) 



and its c.c. 



(5) 

where ci = (-62/2 + Hol^jb^T^ C2 = {1/2Q - b^lDh^T^ 
and C3 = (-&2/2?o - &3/2/5^ + l/2llll)biT^ are complex 
numbers. Note that ci — c\. 

We solve for the Fourier transforms u{k,t) and v{k,t) 
of the fluctuations u(r, t) and u*{r, t), respectively, which, 
from Q satisfy: 



du{k, t) 

dt 
dv(k,t) 

dt 



i^k^ + Ci)u{k) + C2v{k) 

(Ik^ + dHk) + cluik) 



c(fc) 

c*(fc),(6) 



where c(A:) = c{~k) — — Fca X]-exp(«fc-r). Note also that 
u*{—k, t) = v{k, t). The non-zero modes are governed by 



follows similar behavior (Cf. 30]). 

Next, we want to study the effect of spatial expansion. 
We work in new coordinates of rapidity 77 = l/21n(f -|- 
z)/(t — z) and proper time r : 



t — Tcosh{ri) 



Tsinh{ri). 



(11) 



Here we assume the z-axes to be the axes of expansion 
(collision). The change in the corresponding part of the 
Minkowski metric is dt^ — dz"^ = dr^ — r^dr]^, with the 
transverse coordinates left unchanged. Keeping the ra- 
pidity constant results in a constant rate of expansion in 
the system, the speed of expansion being proportional to 
the distance from the collision center: v = z/t, with t the 
time after collision. In the new coordinates the metric is 
similar to the Robertson- Walker metric and is defined as 
ds^ — g^vdx^dxi, ~ dr^ — a^{T)d{Tor])'^ — dx^, i.e. 



g^^ =diag (1, -1, -1, -a^ir)) 



(12) 



where air) — t/tq is the 'scale factor', and tq is the 
parameter controlling the rate of expansion. 

The obvious naive generalization of the model above is 
to substitute 

gijdidj 

in the action , and to consider dynamics in the proper 
frame. 

Eq. ((2Jl now gives 



OT 



-Vll-a'^{T)Vjl + 



dV{l) 



dl 



and its c.c. , 
(13) 
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where V^^ = 9^ + 9^ and V| = 9^. Going through the 
previous development amounts to the naive substitution 
P ^ fc^ + a'^{T)kf, in the eigenvalues: 



This now gives a growth of fluctuations governed by a 
3rd order polynomial in r in the exponent: 



^i,2(fc) = -r 



-{kl + a^T)kl) + c,±\c2\ 



(14) 



u(fc, t) 
v{k, t) 



CiSiexp wi(fc,T')dT' + C2S'2exp / tj2(fc, r')(iT'. 



(15) 



r 



Thus we observe that if 

i(fci+a2(r)fc|)<ci-|c2|, 

there is an explosive growth, much faster than in the 
non-expanding case. Let us for simplicity consider only 
longitudinal modes. We see that the critical mode is 
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a2(r) 



(ci-|c2|). 



(16) 



This is rather remarkable dynamics where the critical 
mode (and all modes at momentum scales below it) is 
moving in time. During initial time there are more scales 
involved but as time progresses only very large scale re- 
gions participate. The prediction of growth under ex- 
pansion by a higher than linear power exponent is tested 
against the simulation results in the following section. 



III. SIMULATION STUDY OF 
DECONFINEMENT DYNAMICS IN SU{3) LGT 



Glauber universality class. The time flow is proportional 
to the link-based heatbath with links visited in random 
order. Therefore, the peak values of the non expan ding 
quench are slightly different from previous studies |3(]l | 
where a heat bath update is performed on links visited 
in systematic order. 

We use a space-time anisotropic lattice in order to be 
able to independently control temperature and expan- 
sion. The anisotropic action is 



S = -/Jj/S^Re 



i>j 



TtU. 



(17) 

where i, j runs over space-like directions, ^ = a/a-r is the 
space-time anisotropy, — G/g^, and = go- ■ gr is the 
anisotropic coupling. Then, by varying two parameters 
it is possible to carry out expansion of the system while 
maintaining constant temperature, or also let the tem- 
perature drop thus allowing for cooling of the expanding 
plasma. 



In this section we present a numerical study of the ef- 
fects of spatial expansion and varying temperature on 
the dynamical evolution following a sudden quench into 
the deconfined phase of the SU{3) gauge theory on the 
lattice. We thus try to mimic conditions encountered in 
heavy-ion collisions. We use periodic boundary condi- 
tions, as we do not study the effects of finite size per 
se. A study of the role of the finiteness of the system 
would certainly be of interest, but is left for a future 
study. Also, in the customary heavy-ion collision picture 
the initial expansion is one-dimensional, becoming three- 
dimensional at later stages of the evolution [s^- Here, as 
we discuss further below, we simplify matters by consid- 
ering a uniform expansion in all spatial directions. This 
is related to isotropic expansion in the proper frame at 
a fixed rapidity value. We first study the effect of such 
expansion; then we add the effect of varying temperature. 

We use a field based heat bath update of the SU{3) 
fields. To accelerate the dynamics there are 2 attempts 
to update the field per sweep. This is somewhat different 
from standard link-based updates; however, it is still in 



A. Spatial expansion 

Hubble-like uniform expansion of the metric amounts 
to varying the space-like lattice spacing a as 



ao(l 



T , 
To' 



(18) 



where r is the proper time variable, and tq is the param- 
eter which controls the rate of expansion. 

We focus here solely on the role of the expansion. So, 
after the initial quench, we keep the temperature T = 
l/lNraT-) constant throughout the expansion by fixing 
the time-like spacing Or- This implies that the anisotropy 
^ = a/ttr traces the changes in a. Thus, assuming zero 
time anisotropy to be f (r = 0) = 1, one has 



^(t)-1 + -. 

To 



(19) 



Next consider the dependence of the space-like lattice 
spacing on the coupling and anisotropy. The one-loop 
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TABLE L The lattice scale parameter A dependence on the 
anisotropy. 





1 


2 


3 


4 


5 


6 


7 




1.000 


0.837 


0.801 


0.798 


0.804 


0.812 


0.820 




8 


9 


10 


11 


12 


13 


14 




0.827 


0.833 


0.838 


0.843 


0.847 


0.851 


0.854 




15 


16 


17 


18 


19 


20 


21 


A(e)/AB 


0.858 


0.860 


0.863 


0.865 


0.867 


0.869 


0.871 



order renormalization group relationship is 

aAiO = cxp{-l/(26o5|)}, (20) 

where bo = 11- S/i^Sn^), and A(^) is dependent on f [iof 
through 

A{0/Ae = exp{-(c,(e) + c,(e))/46o, (21) 

where c^{^) and Ct-(^) are known functions. Inclusion of 
the next order terms adds minor corrections for the range 
of the couplings and anisotropies used and is straightfor- 
ward. It does, however, complicate numerical treatment, 
since it requires a numerical solving of the corresponding 
equation. 

Following the procedure in ^3 compute A(^) for a 
range of anisotropy values as listed in tableU This allows 
us to estimate the necessary time evolution of 



/34(t) =/3(0)-6-26olog 



(1 



To' Ae 



(22) 



It is important to indicate here that a non-perturbative 
further correction 44] need to be applied, when appro- 
priate (see subsection 3.3 below). 



200 
180 
160 
140 
120 
100 
80 
60 
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no expansion 
To=1000 
1(5=500 



not varied; the anisotropy is varied by decreasing the 
time-like lattice spacing. This procedure keeps the cou- 
pling within the scaling window provided the initial cou- 
pling is close to the continuum limit. Here, on the con- 
trary, we keep the time- like spacing constant (or, later, 
slowly increase it) as we vary the space-like coupling. 
This induces changes in the coupling that may drive its 
value out of the scaling regime (/3 ~ 5). Therefore, our 
expansion has to be truncated whenever the value of 
/3 falls below this cut-off value. This condition implies 
that, in order to follow the system evolution for longer 
time, one needs lattices with larger Nr, thus rendering 
the problem more computationally intensive. 

We start simulations on smaller lattices where it is eas- 
ier to gather satisfactory statistics for highly fluctuating 
quantities, such as the structure function. The quench is 
performed from = 5.5 to (3^ = 5.92 on 16'^x4 lattice. 
The latter corresponds to a temperature after the quench 
Tfinai = l.STTc- The phase transition on this lattice at 
^ = 1 is at /? = 5.6902(2). We use jack knife average over 
10 bins, each of 50 configurations. The system is allowed 
to equilibrate for 200 lattice sweep, and then, after per- 
forming a quench, we allow the system to evolve for 800 
sweeps, while measuring several lower modes of the struc- 
ture function. We present here only averages of on-axis 
modes, such as permutations of {n, 0, 0) - this is the n-th 
mode in our notation. The first modes are presented in 
Fig. n We see that expansion significantly enhances the 
response. The faster the expansion rate, the higher are 
the peaks. On the other hand the shift in the location 
of the peaks is not as pronounced, an important point to 
which we return below. 




100 200 300 400 500 600 700 800 



100 200 300 400 500 600 700 800 



FIG. 2: Higher (fe and fcs) structure function modes for non- 
expanding and To = 1000 expanding system. 



FIG. 1: First mode of the structure function at different in- 
verse expansion rates to = 500, 1000, cxd. Quench /Jj — 5.5 
5.92 on 16^x4 lattice. 

Note that in the standard application of anisotropic 
lattices (such as in ji^l) the space- like lattice spacing is 



Next, we look at higher modes of the structure function 
in the cases of no expansion, and expansion at tq = 1000 
- see Fig. |21 We see that the second mode shows behavior 
similar to that observed for the first mode. The differ- 
ence, however, is not that significant. The third mode in 
the expanding system outgrows the corresponding mode 
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in the non-expanding case at early times, but then de- 
creases faster. This is an indication of the shift of the 
critical mode with time, observed in the linear effective 
model (section 2). At later times the transition proceeds 
only through the lower modes. The error bars for some 
of the data points are also presented in Fig. n]and[21 

To make a comparison to the linear effective theory of 
Section 2, we make fits to the exponent for the struc- 
ture function data (Fig. We use a 32'^x4 lattice since 
we know from previous studies [sol Is^ that the linear 
response behavior (pertinent to early times) manifests 
itself better on the larger lattices. Contrary to our ex- 
pectations from the linear theory of section 2, however, 
we find there no substantial change in the behavior of 
the exponent between the expanding and non-expanding 
systems. A fit to 

S{t) - exp(C • r") 

gives a = 1.22 for the expanding system, whereas it gives 
a = 1.18 for the non-expanding system. In Fig. Elwe also 
show fits to S{t) ~ exp(C ■ r), which in fact provides the 
best fit per parameter degree of freedom. Both fits are 
over the t range from to 250. In any case, there is no 
substantial deviation from exponential growth with lin- 
ear T dependence in the exponent. Furthermore, from the 
figures we observe divergence from exponential behavior 
at later times r > 200. Also notable is the enhancement 
of the signal (by a factor ~ 7) in the expanding system as 
compared to the non-expanding system. All this provides 
a manifestation of the limitations of the linear response 
effective models in Section 2. We note that there are 10 
times as many points for the fitting as on the plot. 
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exp(T) 
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50 



100 



150 



200 



250 




50 100 150 200 250 



FIG. 3: Fits to the exponent. Quenches on the 32^x4 lattice 
in non-expanding (top) and tq = 500 expanding (bottom) 
systems. 



B. Expansion accompanied by temperature fall-off 

In this second part of the numerical study we try to 
mimic conditions similar to those in heavy ion collisions 
at RHIC. We want to follow the evolution of the expand- 
ing and cooling QGP. We let the temperature drop as 



(23) 



(1 + T/r^)« 



In the fire-timncl (in the proper frame) the temperature is 
expected to decrease with the proper time r as T ~ -;--i/3 
or slower Q|. We adopt this value of a in the following. 

We first work on smaller lattices for elucidating the 
main physical features. Therefore, we do not set the 
freeze-out condition (see below) and simply take Tq — tq. 
With the a = 1/3 choice of temperature evolution, the 
anisotropy evolves as 



ar) = {i + - 



2/3 



(24) 



and f3 evolves as in (|22ll . We illustrate the effect of tem- 
perature drop on the structure function in Fig. ^ There 



are two basic conclusions suggested by these plots. 

The first is that the temperature evolution drives the 
system back towards the confined phase, but the expan- 
sion tends to prevent it as evidenced by the different peak 
heights. These are then two competing effects that tend 
to cancel, so that the return to structure function equili- 
bration occurs at about the same time (here after about 
700 sweeps) as for the system in the absence of expan- 
sion and temperature fall-off. For the slower expansion 
rates {tq > 1000) and faster temperature fall-offs this 
cancellation effect is even more pronounced. 

The second conclusion is suggested by the fact, also 
present and remarked upon in Fig. ^ that the location 
of the peaks is little affected by the presence of expansion 
and/or temperature evolution. This implies that the sys- 
tem's response to the sudden quench is set by an internal 
dynamics scale that is faster than that of the expansion 
and accompanying temperature fall-off rates considered 
here. After the structure function is past its peak, the 
system is isotropized, and, after a relatively short time, 
any memory of the initial fast, spinodal-like, long range 
response to the violent quench across the deconfinemnt 
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80 




100 200 300 400 500 600 700 800 



FIG. 4: Comparison of first mode of structure functions in- 
cluding drop in temperature on a 16^ x 4 lattice (same quench 
as in Fig. 1). 

transition boundary disappears. The subsequent evolu- 
tion is that of (quasi)equilibrium evolution of the system 
as it expands and cools towards its return to the confine- 
ment phase. 

To elucidate this further we plot the time profile of 
the Polyakov loop average in Fig. \^ for the expanding 
and non-expanding systems as well as with and with- 
out drop in temperature. We now have to use a big- 
ger, 32^ X 8 lattice in order to follow the evolution over 
a longer time interval. First note that pure expansion 

0.6 I 1 , , , 1 1 

expansion and T drop 

expansion, T=const 

0.5 - no expansion, T=const " 

0.4 - 

0.3 - 

0.2 - 



0.1 - 

' 




On the other hand, decreasing temperature counter- 
acts the expansion effect. This is clearly seen in Fig. 
O Eventually, under the combined effects of expansion 
and temperature fall-off, the Polyakov loop expectations 
drops to zero signaling the return of the system to the 
confinement phase. These qualitative features of Fig. [51 
are rather generic, being stable under changes in the ex- 
pansion and temperature fall-off rates (cf. Fig. below). 

The crucial feature characterizing the system's overall 
evolution following the rapid quench into the deconfine- 
ment region is clearly revealed by examining Fig. |S1 in 
conjunction with the plots of the structure function in 
Fig. 2] It is the fact that there are two scales involved in 
this evolution. The first is the scale set by the location of 
the peak of the structure function; the second is the scale 
set by the interval to return to the confinement phase. 
Furthermore, there is wide separation between these two 
scales. As noted above, the location of the peak (~ 300 
in Fig. ^ is very little affected by the conditions of 
expansion and temperature fall-off. It reflects strongly 
coupled dynamics at short time scales driving the expo- 
nential growth of very long range modes leading to 
'isotropization', by which we mean nothing more specific 
than that the full development of these modes appear 
to completely wipe out any remnants of the quenching 
event. Rapid equilibration follows within an interval of 
a few hundreds sweeps (~ 300 — 400 in fig. ^J. The sys- 
tem then continues to evolve in quasi-equilibrium over a 
much longer period (typically of thousands of sweeps as 
in Fig. [S] and Fig. |^ below) expanding and cooling till 
it returns to the confinement phase. This separation of 
time scales is a very general feature over a wide range of 
parameters. 

Having reached this qualitative physical picture, which 
accords well with that experimentally observed, it is in- 
teresting to explore whether it is possible to make an 
estimate of this separation in physical units, and estab- 
lish some correspondence with heavy-ion collision phe- 
nomenology. This we do in the following subsection. 



C. Isotropization-thermalization and chemical 
freeze-out times 



-0.1 I ' ' ' ' ' 1 

500 1000 1500 2000 2500 3000 

FIG. 5: The Polyakov loop evolution after a quench in de- 
creasing temperature on 32^ x 8 lattice (same initial and final 
temperature quench as in Fig. 1). 

drives the Polyakov loop towards larger values, i.e. more 
pronounced deconfined behavior. To get some insight 
into this behavior observe that, as it is evident from (|17|l . 
spatial expansion enhances the timelike part of the action 
(2nd term in the square brackets on the r.h.s.) while sup- 
pressing the spacelike parts ^ij. This tends to increase 
the expectation of timelike Polyakov lines ■ 



The expansion and temperature variation parameters 
(tq, To, Tfinai) determine the precise time interval before 
returning to the confinement phase. Note that our re- 
quirement that the evolution remain inside the scaling 
window limits the time of observation. To deal with this 
and follow the evolution for substantial time we now work 
on a larger, 32'^x8 lattice throughout, and attempt to set 
the parameters so as to reflect conditions in heavy ion 
collision at RHIC. There is strong evidence that the mat- 
ter in the firetunnel reaches temperatures above 2Tc 0. 
For our numerical simulation we quench to Tfinai ~ 3Tc. 
(Note that Tc here means critical temperature of pure 
gluodynamics). The system is equilibrated in the confine- 
ment phase at the same temperature as before « O.STc. 
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For Nt = 8 and ^ = 1 we recalculate the values of corre- 
sponding betas using the two-loop order formula connect- 
ing lattice spacing and the coupling, and reweigh it with 
non-perturbative correction factor above Tc |4J| . Such a 
correction is needed here since the bare gauge coupling is 
typically of order one for our range of betas. The value of 
f3c is known from the lattice Polyakov loop susceptibility 
study ,4^. We get 



initial =5.90 O.reTe 

Pc = 6.0625 Tc 
/3fi„ai = 6.85 2.95r, , 



(25) 



where 'initial' and 
quench. 



■final' refer to before and after the 
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FIG. 6; The Polyakov loop evolution after a quench to 2.95Tc, 
TO = 1000 

At high enough temperatures the system is close to an 
ideal gas of quasiparticles. Lattice data suggest that the 
equation of state does not actually deviate much from 
that of the ideal gas for all temperatures down to Tc, 
and is conveniently modeled as such in hydrodynamic de- 
scriptions. Hence, the energy and entropy densities scale 
as e ~ T'* and s ^ T^, respectively. Assuming adiabatic 
expansion the entropy per unit of rapidity is conserved. 
In real fire-tunnel evolution the initial expansion is one 
dimensional [39l |. switching to three dimensional expan- 
sion at later times. This corresponds to T ~ r~^/^ at 
the early and mid time and T ^ t^^ at later time. The 
T ~ T~^^^ behavior is in fact seen in the hydrodynamic 
evolution almost down to Tc jlf, ^ 3] j a-nd is the only one 
considered in the following. We are thus led to model- 
ing of the spatial expansion and temperature drop on the 
lattice in terms of the two parameters tq and as: 



ao(l 
ao(l 



To 

^11/3 



ao(l + 2/-)l/^ 

To 



(26) 



(27) 



where y = tq/tq is the ratio of the 'speeds'. The evolution 
of the anisotropy then is 



e(r) 



1 + f 



(1 + 2/^)^/^ 



(28) 



As the chemical freeze-out temperature we choose 
Tfo — Tc 0, 0, Before this freeze-out the plasma 
undergoes an a;-fold expansion 



ao 



It. 

To 



On the other hand 



Tfo 



Tfrnal (l+Z/^)l/3 

From these two we get 

{Tfinal/TfoY - 1 



y = 



X - 1 



(29) 



(30) 



(31) 



Hydrodynamical model phenomenology yields all re- 
quired parameters 0- Thus we have a: ~ 9 before the 
conversion to the confined phase is completed, and the 
freeze-out is observed. This corresponds to y ~ 3.25. 

In Fig. Elwe plot the time profile of the Polyakov loop 
average in expanding and non-expanding systems with 
and without drop in temperature with these values of x, 
y, and tq = 1000. One again observes the same features 
discussed above in connection with Fig. |31 




FIG. 7: The first mode of structure function 
expanding and constant T system after a quench. 



for 



In Fig. [7|we show the lowest mode structure function 
in the case of no expansion and temperature variation. 
The expanding-variable temperature case is very similar. 
We see that at r = 400 sweeps the system has reached 
the peak. This is the point of isotropisation of the system 
when the long range fluctuations reach through the sys- 
tem. It is again to be contrasted to the much longer times 
needed to return to confinement (vanishing Polyakov loop 
expectation) . 
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FIG. 8: The Polyakov loop evolution after a quench to 2.95Tc 
for various values of tq and x = 9 (top panel), a; = 10 (bottom 
panel) . 



To explore this difference in detail, in Fig. |51 we 
present the Polyakov loop average evolution at different 
expansion rates. The figure on the top panel uses again 
X — 9, so the chemical freeze-out point on this plot is 
at t/tq = (a; — 1) = 8. We notice that in the range of 
To ~ 400 — 600 the Polyakov loop expectation value gets 
close to zero just before this point. Actually, this should 
be considered as a lower bound on the range of tq's since 
this is a first order transition (or a rapid crossover in the 
presence of fermions). Thus, there is a latent heat pe- 
riod during which the system lingers, with Polyakov loop 
expectations close to zero, before it is fully converted 
to the confined phase. The values tq = 300 and less 
(higher 'speeds') do not show this behavior and there- 
fore do not lead to a consistent picture. Taking then 
the range of tq > 400 — 600 gives an interval to freeze- 
out of > 3200 — 4800 sweeps, which in turn corresponds, 
from phenomenology, to ~ 9/m/c. This allows an up- 
per bound estimate of the time around the peak of the 
structure fmiction in physical units: < 0.75 — 1.1 fm/c. 

In the plot on the bottom panel in Fig. |Slwe use the 
same quench to 2.95Tc but now take x = 10, which gives 



y = 2.89. This corresponds to a somewhat longer life- 
time for the deconfined plasma. One sees that a value of 
To ^ 700 — 800 gives a good fit to complete conversion 
to the confined phase by the time the freeze-out point is 
reached. This in turn gives < 0.56 — 0.63 fm/c as an 
upper bound estimate for the time around the peak of 
the structure function. There is of course an inherent 
uncertainty here as to what to take for the appropriate 
phenomenological value for the lifetime in physical units 
since we do not include fermions in out simulations; and 
fermions become important at the late stage of conver- 
sion back to confinement. In Fig. |Slwc follow the curves 
as far as possible before running out of the scaling regime 
as explained in Section 3.1 above. Longer lifetimes result 
into even shorter isotropization times. To consider some- 
what larger x values, however, we need larger lattices 
than those employed in this study. But the main mes- 
sage extracted from the present decimations should be 
clear. The robust separation between the 'fast' dynam- 
ics scale of the exponential growth of the spinodal-like 
response to the sudden quench and that set by return 
to confinement is such that isotropization times well in- 
side 1 fm/c result naturally for any reasonable choice of 
parameters. 



IV. CONCLUSIONS 

In this paper we studied the response of the pure SU (3) 
gauge theory to a rapid quench from its confined to its de- 
confined phase. In a series of simulations wc followed the 
subsequent evolution of the system under varying con- 
ditions of temperature fall-off and/or spatial expansion. 
These conditions were chosen so as to reflect the type of 
variations presumed to hold in heavy ion collisions. Our 
main finding is that there are two distinct scales char- 
acterizing this evolution. There is one scale set by the 
development of very long range modes continuously from 
zero to their maximum amplitude over a short time in- 
terval. These modes, manifested in the response of the 
structure functions to the quench, drive isotropization 
and return to thermalization shortly afterwards. The 
scale set by this 'fast' dynamics is little affected by condi- 
tions of expansion and temperature variations. The sec- 
ond scale is set by the time interval to return to the con- 
finement phase under quasi-equilibrium evolution, and 
is affected by the details of expansion and temperature 
fall-off. There is a robust wide separation between the 
two scales, which, translated to physical units under rea- 
sonable assumptions about the lifetime of the plasma, 
gives estimates of the 'fast' dynamics in the range of 
0.5 - 1 fm/c. 

There is a number of directions in which this study can 
be extended. The formalism used above can be extended 
to treat also anisotropy among the different space direc- 
tions, thus allowing a more 'realistic' treatment of spatial 
expansion. Whether this makes much of difference, how- 
ever, remains to be seen (cf. footnote j42i|). Another 
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interesting question is that concerning the effect of hav- 
ing a really finite physical system. In our simulations we, 
as usual, employed periodic boundary conditions. But 
one may explore different ones that would mimic a finite 
rather than an infinite system. The other obvious ex- 
tension is the inclusion of fermions. The deconfinement 
transition is driven by gluonic dynamics. Fermions, how- 
ever, are expected to contribute more substantially at the 
late stage of return to confinement and hadronization. 
Thus our study applies strictly to the gluonic plasma. 
Nonetheless, the qualitative picture of the separation of 
scales found above should not be affected in an essential 
way by the inclusions of fermions, though quantitative 



details related to the exact lifetime of the plasma, etc 
certainly will. 
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